/* ****************************************************************** **
**    OpenSees - Open System for Earthquake Engineering Simulation    **
**          Pacific Earthquake Engineering Research Center            **
**                                                                    **
**                                                                    **
** (C) Copyright 1999, The Regents of the University of California    **
** All Rights Reserved.                                               **
**                                                                    **
** Commercial use of this program without express permission of the   **
** University of California, Berkeley, is strictly prohibited.  See   **
** file 'COPYRIGHT'  in main directory for information on usage and   **
** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
**                                                                    **
** Developed by:                                                      **
**   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
**   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
**   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
**                                                                    **
** ****************************************************************** */
                                                                        
// $Revision: 1.10 $
// $Date: 2007-04-02 23:41:13 $
// $Source: /usr/local/cvs/OpenSees/SRC/analysis/algorithm/equiSolnAlgo/NewtonRaphson.cpp,v $
                                                                        
                                                                        
// File: ~/OOP/analysis/algorithm/NewtonRaphson.C 
// 
// Written: fmk 
// Created: Sun Sept 15 15:06:47: 1996 
// Revision: A 
//

// Description: This file contains the class definition for 
// NewtonRaphson. NewtonRaphson is a class which uses the
// Newton-Raphson solution algorithm
// to solve the equations. No member functions are declared as virtual as 
// it is not expected that this class will be subclassed.
// 
// What: "@(#)NewtonRaphson.C, revA"

#include <NewtonRaphson.h>
#include <AnalysisModel.h>
#include <StaticAnalysis.h>
#include <IncrementalIntegrator.h>
#include <LinearSOE.h>
#include <Channel.h>
#include <FEM_ObjectBroker.h>
#include <ConvergenceTest.h>
#include <ID.h>
#include <elementAPI.h>
#include <string>

// Forward declaration for GPUSolver
class GPUSolver;

#include <chrono>
using namespace std::chrono;

	// Global timing output control
bool enableTimingOutput = 0;  // ← 设置为0关闭输出

// ====== Performance Optimization: CPU formTangent skip counter ======
static int cpuFormTangentCalls = 0;
static int gpuFormTangentSkips = 0;
// ====== End of Performance Optimization ======

void* OPS_NewtonRaphsonAlgorithm()
{
    int formTangent = CURRENT_TANGENT;
    double iFactor = 0;
    double cFactor = 1;

    while (OPS_GetNumRemainingInputArgs() > 0) {
      const char* type = OPS_GetString();
      if(strcmp(type,"-secant")==0 || strcmp(type,"-Secant")==0) {
	formTangent = CURRENT_SECANT;
	iFactor = 0;
	cFactor = 1.0;
      } else if(strcmp(type,"-initial")==0 || strcmp(type,"-Initial")==0) {
	formTangent = INITIAL_TANGENT;
	iFactor = 1.;
	cFactor = 0;
      } else if(strcmp(type,"-intialThenCurrent")==0 || strcmp(type,"-intialCurrent")==0) {
	formTangent = INITIAL_THEN_CURRENT_TANGENT;
	iFactor = 0;
	cFactor = 1.0;
      } else if(strcmp(type,"-hall")==0 || strcmp(type,"-Hall")==0) {
	formTangent = HALL_TANGENT;
	iFactor = 0.1;
	cFactor = 0.9;
	if (OPS_GetNumRemainingInputArgs() == 2) {
	  double data[2];
	  int numData = 2;
	  if(OPS_GetDoubleInput(&numData,&data[0]) < 0) {
	    opserr << "WARNING invalid data reading 2 hall factors\n";
	    return 0;
	  }
	  iFactor = data[0];
	  cFactor = data[1];
	}
      }
    }

    return new NewtonRaphson(formTangent, iFactor, cFactor);

}

// Constructor
NewtonRaphson::NewtonRaphson(int theTangentToUse, double iFact, double cFact)
:EquiSolnAlgo(EquiALGORITHM_TAGS_NewtonRaphson),
 tangent(theTangentToUse), iFactor(iFact), cFactor(cFact)
{

}

NewtonRaphson::NewtonRaphson()
	:EquiSolnAlgo(EquiALGORITHM_TAGS_NewtonRaphson),
	tangent(CURRENT_TANGENT), iFactor(0.), cFactor(1.)
{

}


NewtonRaphson::NewtonRaphson(ConvergenceTest &theT, int theTangentToUse, double iFact, double cFact)
:EquiSolnAlgo(EquiALGORITHM_TAGS_NewtonRaphson),
 tangent(theTangentToUse), iFactor(iFact), cFactor(cFact)
{

}

// Destructor
NewtonRaphson::~NewtonRaphson()
{
  

}


int NewtonRaphson::solveCurrentStep(void)
{
	auto t_start = high_resolution_clock::now();

	AnalysisModel* theAnaModel = this->getAnalysisModelPtr();
	IncrementalIntegrator* theIntegrator = this->getIncrementalIntegratorPtr();
	LinearSOE* theSOE = this->getLinearSOEptr();

	if ((theAnaModel == 0) || (theIntegrator == 0) || (theSOE == 0)
		|| (theTest == 0)) {
		opserr << "WARNING NewtonRaphson::solveCurrentStep() - setLinks() has";
		opserr << " not been called - or no ConvergenceTest has been set\n";
		return -5;
	}

	auto t1 = high_resolution_clock::now();
	if (theIntegrator->formUnbalance() < 0) {
		opserr << "WARNING NewtonRaphson::solveCurrentStep() -";
		opserr << "the Integrator failed in formUnbalance()\n";
		return -2;
	}
	auto t2 = high_resolution_clock::now();

	theTest->setEquiSolnAlgo(*this);
	if (theTest->start() < 0) {
		opserr << "NewtonRaphson::solveCurrentStep() -";
		opserr << "the ConvergenceTest object failed in start()\n";
		return -3;
	}

	int result = -1;
	numIterations = 0;

	// ͳ���ܺ�ʱ
	long long total_formTangent = 0;
	long long total_solve = 0;
	long long total_update = 0;
	long long total_formUnbalance = 0;
	long long total_test = 0;

	do {
		auto t_ft1 = high_resolution_clock::now();

		// ====== Performance Optimization: Skip CPU formTangent for GPU Solver ======
		// 检测是否是GPU求解器，如果启用了GPU组装，则跳过CPU端formTangent
		// 使用类名判断而非动态转换，避免头文件依赖
		bool skipFormTangent = false;
		const char* soeName = theSOE->getClassType();
		if (soeName != nullptr && strcmp(soeName, "GPUSolver") == 0) {
			skipFormTangent = true;
			gpuFormTangentSkips++;
			if (gpuFormTangentSkips == 1) {
				printf("[CPU_OPTIMIZATION] GPU Solver detected - skipping CPU formTangent\n");
			}
		}

		if (!skipFormTangent) {
			// 传统求解器：执行CPU端formTangent
			cpuFormTangentCalls++;
			// --- formTangent ---
			if (tangent == INITIAL_THEN_CURRENT_TANGENT) {
				if (numIterations == 0) {
					SOLUTION_ALGORITHM_tangentFlag = INITIAL_TANGENT;
					if (theIntegrator->formTangent(INITIAL_TANGENT) < 0) {
						opserr << "WARNING NewtonRaphson::solveCurrentStep() -";
						opserr << "the Integrator failed in formTangent()\n";
						return -1;
					}
				}
				else {
					SOLUTION_ALGORITHM_tangentFlag = CURRENT_TANGENT;
					if (theIntegrator->formTangent(CURRENT_TANGENT) < 0) {
						opserr << "WARNING NewtonRaphson::solveCurrentStep() -";
						opserr << "the Integrator failed in formTangent()\n";
						return -1;
					}
				}
			}
			else {
				SOLUTION_ALGORITHM_tangentFlag = tangent;
				if (theIntegrator->formTangent(tangent, iFactor, cFactor) < 0) {
					opserr << "WARNING NewtonRaphson::solveCurrentStep() -";
					opserr << "the Integrator failed in formTangent()\n";
					return -1;
				}
			}
		}
		else {
			// GPU求解器：跳过CPU端formTangent，仅设置tangent flag（供GPU solver使用）
			if (tangent == INITIAL_THEN_CURRENT_TANGENT) {
				SOLUTION_ALGORITHM_tangentFlag = (numIterations == 0) ? INITIAL_TANGENT : CURRENT_TANGENT;
			}
			else {
				SOLUTION_ALGORITHM_tangentFlag = tangent;
			}
		}
		// ====== End of Performance Optimization ======

		auto t_ft2 = high_resolution_clock::now();

		// --- solve ---
		if (theSOE->solve() < 0) {
			opserr << "WARNING NewtonRaphson::solveCurrentStep() -";
			opserr << "the LinearSysOfEqn failed in solve()\n";
			return -3;
		}
		auto t_solve = high_resolution_clock::now();

		// --- update ---
		if (theIntegrator->update(theSOE->getX()) < 0) {
			opserr << "WARNING NewtonRaphson::solveCurrentStep() -";
			opserr << "the Integrator failed in update()\n";
			return -4;
		}
		auto t_update = high_resolution_clock::now();

		// --- formUnbalance ---
		if (theIntegrator->formUnbalance() < 0) {
			opserr << "WARNING NewtonRaphson::solveCurrentStep() -";
			opserr << "the Integrator failed in formUnbalance()\n";
			return -2;
		}
		auto t_fu = high_resolution_clock::now();

		// --- test ---
		result = theTest->test();
		auto t_test = high_resolution_clock::now();

		numIterations++;
		this->record(numIterations);

		// �ۼ�ÿ�ε������׶κ�ʱ
		total_formTangent += duration_cast<milliseconds>(t_ft2 - t_ft1).count();
		total_solve += duration_cast<milliseconds>(t_solve - t_ft2).count();
		total_update += duration_cast<milliseconds>(t_update - t_solve).count();
		total_formUnbalance += duration_cast<milliseconds>(t_fu - t_update).count();
		total_test += duration_cast<milliseconds>(t_test - t_fu).count();

	} while (result == -1);

	if (result == -2) {
		opserr << "NewtonRaphson::solveCurrentStep() -";
		opserr << "the ConvergenceTest object failed in test()\n";
		return -3;
	}

	auto t_end = high_resolution_clock::now();

	// 性能总耗时统计
	if (enableTimingOutput) {
		printf("formUnbalance:  %lld ms\n", duration_cast<milliseconds>(t2 - t1).count() + total_formUnbalance);
		printf("formTangent:    %lld ms (CPU calls: %d, GPU skips: %d)\n",
			   total_formTangent, cpuFormTangentCalls, gpuFormTangentSkips);
		printf("solve:          %lld ms\n", total_solve);
		printf("update:         %lld ms\n", total_update);
		printf("test:           %lld ms\n", total_test);
		printf("total:          %lld ms\n", duration_cast<milliseconds>(t_end - t_start).count());

		// 重置计数器（每个step后）
		cpuFormTangentCalls = 0;
		gpuFormTangentSkips = 0;
	}

	return result;
}


int
NewtonRaphson::sendSelf(int cTag, Channel &theChannel)
{
  static Vector data(3);
  data(0) = tangent;
  data(1) = iFactor;
  data(2) = cFactor;
  return theChannel.sendVector(this->getDbTag(), cTag, data);
}

int
NewtonRaphson::recvSelf(int cTag, 
			Channel &theChannel, 
			FEM_ObjectBroker &theBroker)
{
  static Vector data(3);
  theChannel.recvVector(this->getDbTag(), cTag, data);
  tangent = int(data(0));
  iFactor = data(1);
  cFactor = data(2);
  return 0;
}


void
NewtonRaphson::Print(OPS_Stream &s, int flag)
{
  if (flag == 0) {
    s << "NewtonRaphson" << endln;
  }
}


int
NewtonRaphson::getNumIterations(void)
{
  return numIterations;
}


